function [w]=funcEM(const,h,tmax,winit)
t_i = length(0:h:tmax);
w = zeros(3,t_i);   % Preallocate w
w_= zeros(3,t_i);   % Preallocate w_
w(:,1) = winit; % Initial values
for i = 2:t_i
    w_(:,i) = w(:,i-1)+h*func2U8(w(:,i-1),const); % Predictor
    w(:,i) = w(:,i-1)+(h/2)*(func2U8(w(:,i-1),const)+func2U8(w_(:,i),const)); % Corrector
end